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Starting from an exact lower bound on the imaginary-time propagator, we present a Path- 
Integral Quantum Monte Carlo method that can handle singular attractive potentials. We illustrate 
the basic ideas of this Quantum Monte Carlo algorithm by simulating the ground state of hydrogen 
and helium. 
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I. INTRODUCTION 



Quantum Monte Carlo (QMC) simulation is a powerful method for computing the ground state and non-zero 
■ temperature properties of quantum many-body systems Klpl. There are two fundamental problems that limit the 
application of these methods. The first and most important is the minus-sign problem on which we have nothing to 
O |. say in this paper, see however The second problem arises if one would like to simulate systems with attractive 

singular potentials, the Coulomb interaction being the prime example. The purpose of this paper is to present an 
approach that solves the latter problem, in a form that fits rather naturally in the standard Path Integral QMC 
Q ■ (PIQMC) approach and leaves a lot of room for further systematic improvements. 

Let us first recapitulate the basic steps of the procedure to set up a PIQMC simulation. Writing K and V for the 
kinetic and potential energy respectively the first step is to approximate the imaginary-time propagator by aproduct 
* 55 ' °f short-time imaginary-time propagators. The standard approach is to invoke the Trotter-Suzuki formula (^^] 
>»■ 

O l' m— »oc V / 

to construct a sequence of systematic approximations Z m to the partition function Z 

> 

l/V Z = Tr exp(-(3H) = lim Z m (2) 



o 

O ' J n=l 



/III- 
dm ■ ■■dr m J{{r n \e-e K '™\r n+x )e-P v ^)/™, (3) 



where r m+1 = m and use has been made of the fact that the potential energy is diagonal in the coordinate represen- 
tation. Taking the limit m — > oo, (3) yields the Feynman path integral [0] for a system with Hamiltonian H = K + V. 
Q ' Expression (3) is the starting point for PIQMC simulation. 

In the case the attractive Coulomb interaction, it is easy to see why the standard PIQMC approach fails. Let us 
take the hydrogen atom as an example. The Hamiltonian reads 

B<: H = -^- q -, (4) 

> ; 2M r ' w 

k> , where q denotes the charge of the electron and M = m e /(l + m e /m p ), m e (m p ) being the mass of the electron 
; (proton). Replacing the imaginary-time free-particle propagator in (3) by its explicit, exact expression 



-f>K/ my) = ( -^V /2 exp (_ "M|*-/I a ) , (5) 



we obtain 



mM ^ 3m/2 



Zm = ( „ ^2 ) I dm--- dr m exp 
,27T pn 



mM x - 



2/3h 2 



r n - r n+ i) 2 



71=1 



exp 



n—l 



(6) 



PIQMC calculates ratio of integrals such as (6) by using a Monte Carlo procedure to generate the coordinates 
{m, . . . , r m }. The integrand in (6) serves as the weight for the importance sampling process. As the latter tends to 
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maximize the integrand, it is clear that because of the factors exp (+f3q 2 m~ 1 r~ 1 ), the points {ri, . . . ,r m } will, after 
a few steps, end up very close to the origin. In the case of a singular, attractive potential importance sampling based 
on (6) fails. Using instead of the simplest Trotter-Suzuki formula (1) a more sophisticated one || only makes things 
worse because these hybrid product formulae contain derivatives of the potential with respect to the coordinates. 

The problem encountered in setting up a PIQMC scheme for models with a singular, attractive potential is just a 
signature of the fundamental difficulties that arise when one tries to define the Feynman path integral for the hydrogen 
atom . The formal solution to this problem is known . It is rather complicated and not easy to incorporate 

in a PIQMC simulation. 

In spirit the method proposed in this paper is similar to the one used to solve the hydrogen path integral: Use 
the quantum fluctuations to smear out the singularity of the potential. Mathematically we implement this idea by 
applying Jensen's inequality to the propagator jl2| ]. Applications of the Feynman path- integral formalism are often 
based on a combination of Jensen's inequality and a variational approach |^|,[l0| so it is not a surprise that similar 
tricks may work for PIQMC as well. 

The paper is organized as follows. In Section 2 we give a simple derivation of an exact lower bound on the 
imaginary-time propagator. This inequality naturally defines a sequence of systematic approximations Z m to the 
partition function. Although each Z m looks very similar to Z mi the former can be used for PIQMC with attractive, 
singular potentials. For pedagogical reasons, in Section 3 we illustrate the approach by presenting an analytical 
treatment of the harmonic oscillator. In Section 4 we give the explicit form of the approximate propagator for the 
attractive Coulomb potential and present PIQMC results for the ground state of the hydrogen and helium atom. 

II. LOWER BOUND ON THE PROPAGATOR 

Consider a system with Hamiltonian H = K + V and a complete set of states {|a;)} that diagonalizes the hermitian 
operator V. In the case that V contains a singular attractive part we replace V = lim e ^o U e by a regular V e {x) > — oo 
and take the limit e — > at the end of the calculation. Using the Trotter-Suzuki formula we can write 

(x\e- T ( K+v ^\x') = lim (x\ (e- TK/m e- rV ' /m ) m \x'), (7) 

/m 
dx x ■■■dx n n<^|e- rK/m |^+i>^ ry£( * l)/m , (8) 

i—l 

- lim I 1 a x nlU =1 {x,\e \x l+1 )e dx,---dx Y\ (x\ e - rK/m \x^)(9) 

™-oc jdx 1 ---dx n n? =1 (^- TK/m \^ + i) J dx ™n^ e ^+^> 

If (x\e~ rK \x t ) > for all r, x and x' \ the function 

P({xi}) = Y[(xi\e- TK / m \x i+1 ) / / dx 1 ---dx n ]l(x l \e- TK / m \x l+1 ), (10) 
i=i I •' i=i 

is a proper probability density. Clearly (10) is of the form J dx\ ■ ■ ■ dx n p({xi})f({xi}) so that we can apply Jensen's 
inequality 

dxx--- dxnpHx^e 9 ^^ > exp (J dx x --- dx„p({x. i }) 5 ({x l })^ , (11) 
and obtain 

(x\e-^ K+ ^\x>) > (x\e-^\ x >) lim expf-^f^ f dx x ■ ■ ■ dx m IE=i(*»£*' m \ x » + i) \ (12) 

nwoo \ m ^— ' J {x\e \x ) I 



( T m f 

> (x\e~ TK \x') lim exp >^ / dx, 

m^oo \ m * — ' / 



{x\e- rK \x') 
{x\e TK/m \xi)V e {x i )(x i \ e TK ' m \x') 



(x\e~ TK \x') J ' 

For m — » oo, the sum over n can be replaced by an integral over imaginary time. Finally we let e — > and obtain E 



(13) 
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(x\e- T(K+v) \x') > (x\e- TK \x') 



exp 



du 



(x\e- uK Ve-<- T -^ K \x') 



(x e 



-tKWi\ 



(14) 



Note that l.h.s of (14) reduces to the standard, symmetrized Trotter-Suzuki formula approximation |l^|l4|] if we 
replace the integral over u by a two-point trapezium-rule approximation. This replacement also changes the direction 
of inequality as can been seen directly from the upperbound |12| 



(x\e- T{K+v) \x'} < (x\e- TK \x')exp\- duln 



{x\e- uK e- TV e-^-^ K \x') 

-tK\„i\ 



< {x\e~ T *\x')e 



tK\I\~tV{x) 



Expression (14) can be used to define a new type of approximant to the partition function namely 



(15) 



M 



3m/2 



/m 
dr 1 .. . dr m exp 



M 
2?f? 



'n+1) 



du 



(r n \e- uK Ve-^ K \r n+1 ) 
(r n \e- TK \r n+ i) 



(16) 



where r = /3/m. The simplest approximant Z\ corresponds to the Feynman's variational approximation to the full 
Feynman path integral p|JIO[|. The main difference between (3) and (16) is that the bare potential e~ rV ^ is replaced 
by an effective potential that is obtained by convoluting the bare potential and free-particle propagators e~ uK and 
6 -(t-u)K_ Convolution smears out singularities. As we show below, in the case of the attractive Coulomb interaction 
expression (14) is finite, for any choice of x and x' . For the approximants Z m to be useful in PIQMC, it is necessary 
that the integral over u can be done efficiently. In the next two sections we show how this can be done. 

III. ILLUSTRATIVE EXAMPLE 

It is instructive to have at least one example for which the details can be worked out analytically, without actually 
using PIQMC. Not surprisingly this program can be carried out for the harmonic oscillator. For notational convenience 
we will consider the one-dimensional model Hamiltonian H = K + V, with K = — (h 2 /2M)d 2 /dx 2 and V = Mus 2 x 2 . 
Calculating the matrix element {x\e~ uK Ve~^ T ~ w ^ K \x') in (16) is a straightforward excercise in perfoming Gaussian 
integrals Jig]. We obtain 



Z m — 



mM 
2n(3h 2 



i/2 



dxi 



■ dx m Yl 



exp 



mM _ , 2 (3Mcu 2 2 2 $h 



2p?r 



6m 



2mM' 



.(17) 



The integrand in (17) is a quadratic form and can be diagonalized by a Fourier transformation with respect to the 
index n. Evaluation of the resulting Gaussian integrals yields 



r\ — m/2 

m. — ^ 



exp 



p 2 h 2 Lu 2 ^ m -' 



12m 



n 



3m 



(3 2 h z uo 2 ( 2 n z uj 2 \ /27m 
1 + — 1 - — cos 



6 m 



-1/2 



(18) 



Taking the partial derivative of — In Z m with respect to (3 gives the corresponding approximation to the energy: 



2, ,2 



E m = 



6m 



1 



m— 1 



2 + cos (2im/m) 



~^ 1 — cos (2nn/m) + /3 2 h 2 uj 2 (2 + cos(27rn/m)) /6m 
For comparison, if we use of the standard Trotter-Suzuki formula we obtain @ 

1 



I3h 2 u 2 ^ 



2m 2 ^ 1 - cos (27m/m) + (3 2 h 2 uj 2 /2m 2 



(19) 



(20) 



In Table 1 we present numerical results obtained from (19) and (20) and compare with the exact value of the energy 
E = (h,uj/2)coth((3huj/2)). Note that the average of the two approximations, i.e. (E m + E m )/2, is remarkably close 
to the exact value E, an observation for which we have no mathematical justification at this time. 
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TABLE I. Numerical results for the exact energy of the harmonic oscillator (E), and approximations based on (19) (E m ) 
and (20) (E m ). We use units such that hco = 1 and j3 is dimensionless. 





m 


E m 


E 


E m 


1 


1 


1.00000 


1.08198 


1.16668 




10 


1.08101 


1.08198 


1.08292 




50 


1.08194 


1.08198 


1.08202 




100 


1.08197 


1.08198 


1.08199 




500 


1.08198 


1.08198 


1.08198 


5 


1 


0.20000 


0.50678 


1.03333 




10 


0.49199 


0.50678 


0.51938 




50 


0.50617 


0.50678 


0.50694 




100 


0.50678 


0.50678 


0.50679 




500 


0.50678 


0.50678 


0.50679 


10 


1 


0.10000 


0.50005 


1.76667 




10 


0.44273 


0.50005 


0.54316 




50 


0.49757 


0.50005 


0.50234 




100 


0.49942 


0.50005 


0.50064 




500 


0.50002 


0.50005 


0.50007 
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IV. ATTRACTIVE COULOMB POTENTIAL 



As a second example we will consider a neutral system consisting of two electrons with opposite spin and a nucleus. 
The Hamiltonian reads [fl6| . [i"7[ 

_ H 2 2 h 2 2 q 2 q 2 2q 2 
H --^^-2M 2 V2 -W\~W\ + V^\' (21) 

where the vectors r\ and r 2 describe the position of the two electrons, with the nucleus placed in the origin. It is 
convenient to introduce the notation Ki = — D i Vf, Di = h 2 /2Mi, Vj = V(ri), V12 — V{r\ — r 2 ), and V(r) = q 2 /\r\, 
for i = 1,2. Application of inequality (14) requires the evaluation of 

, 1 \ / T rf U (r 1 r 2 |e-^+^)(y 1 + ^-2^ 12 )e-( T -«)(^+^)|rir 2 } 
I( ri ,r 2 , ri ,r 2 ) = {ri r 2 \e-^^)Kr>) (22) 

Jo duinle-^Vie-^-^lri) ^ du(r 2 \e- uK ^V 2 er^- u ^\r' 2 ) 



(n|e-^i)|ri) (r 2 \e- TK i\r 2 ) 

fl dul'i 

+2- 



JZdu{r 1 r 2 \e-»l Kl+K *)Vi2e-t T -»X K i+ K '^r{r' 2 ) 



where we made use of the fact that [K\, V 2 ] = [K 2 , V\\ — 0. It is sufficient to consider the last term of (22). Inserting 
a complete set of states for both particles we obtain 

f[ du f dr'j f dr^r ir2 \e-^+ K ^\r>>r>>)V(r>> r^{ry<\e-^){K 1+ K 2 ) 
h2(rur 2 , ri ,r 2 ) = (r^e-^+K^) • ( 23 ) 

Inserting the explicit expression for the free-particle propagator (5), a straightforward manipulation of the Gaussian 
integrals in (23) gives 



h2(n,r 2 ,r[y 2 ,D)=J o du j dr ^ <T T _ u)D } V(r)exp{ 



( r™ — (t — u){r\ — r 2 ) — u(r[ — r' 2 )] 2 



rr 



4mt(t — u)D 



, (24) 



where D = D\ + D 2 

In the case of the Coulomb potential, the integral over r can be evaluated analytically by changing to spherical 
coordinates. The remaining integral over u is calculated numerically. In practice, it is expedient to replace the 
integration over u by an integration over an angle. An expression that is adequate for numerical purposes is 

, ^ erf [(4TD)- 1 / 2 1 (ri - rg) tan <j> + (r[ - r> 2 ) cot 0|] 

/i 2 (ri,r 2 ,r 1 ,r 2 ,£)) =2rq / defy — r — — . (25) 

Jo \(ri ~r 2 )tan0+ {r' x -r 2 )cot0| 

It is easy to check that Ix2(fx,T 2 , r[, r' 2 , D) is finite. The expressions for the first and second contributions in (22) can 
be obtained from (25) by putting (D 2 , r 2 , r 2 ) and (Dx, n, r[) equal to zero, i.e. Ix(rx, r[, Dx) = Ixi{tx, 0, r[, 0, Dx) 
zndI 2 {r 2 ,r 2 ,D 2 )=I 12 (0,r 2 ,0,r' 2 ,D 2 ). 

For the helium atom M = Mi = Af 2 , and the m-th approximant to the partition function reads 

3m 



^ = (i^y / ^•••^^•••^ eX p|-^£[( r «- r »+l) 2 + K- r "n+l) 2 ]| ( 26 ) 

xexpjr^ [h(r n ,r n+ i,Di) + I 2 (r' n ,r' n+1 , Dx) - 27i 2 (r n , r n+1 , r' n , r' n+1 , 2D X ) j , 
whereas in the case of the hydrogen atom we have 
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zH 



M s 3m/2 . | m '" I 

J / dri . . . dr m exp < -^-^2 ^(r„ - r„+i) 2 + Ii(r„,r„ + i, £>i) > , (27) 

^ \ 71=1 n—1 ) 



2tttH 2 



with t = (3/m. As the integrands in (26) and (27) are always finite, expressions (26) and (27) can be used perform 
PIQMC simulations. 

In the path integral formalism the ground state energy is obtained by letting (3 — > oo and (3/m — > 0, i.e. E = 
lim^—too lirng/m^o E m . Of course, in numerical work, taking one or both these limits is impossible. In Tables 2 
and 3 we present numerical results of PIQMC estimates of the ground state energy E of the hydrogen and helium 
atom. These results have been obtained from five statistically independent simulations of 100000 Monte Carlo steps 
per degree of freedom each. The systematic errors due to the discretization of the path integral are hidden in the 
statistical noise. The PIQMC procedure we have used is standard [|]J^] except for a trick we have used to improve 
the efficiency of sampling the paths, details of which are given in the appendix. Although a ground state calculation 
pushes the PIQMC method to point of becoming rather inefficient, the numerical results are in satisfactory agreement 
with the known values. 



V. DISCUSSION 



We have show that is possible to perform PIQMC simulations for quantum systems with attractive Coulomb 
potentials. Instead of the conventional Trotter-Suzuki formula approach one can use (16) to construct a path integral 
that is free of singularities. In practice, a numerical calculation of the latter requires only minor modifications of a 
standard PIQMC code. 

The efficiency of the PIQMC method describe above can be improved with relatively modest efforts. Instead of using 
the free-particle propagator K, we are free to pick any other model Hamiltonian Hq for which the matrix elements of 
e -Tffo are p 0S iti ve and integrals involving these matrix elements are known analytically. An obvious choice would be 
to take for Hq a set of harmonic oscillators. The matrix elements of e~ rH ° are Gaussians and hence the conditions 
used to derive (14) are satisfied. If necessary the approximant Z m can be improved further by optimization of the 
parameters of the oscillators. For m = 1 this approach is identical to the variational method proposed by Feynman 
and Kleinert Jl§|-|2~i}| and independently by Giachetti and Tognetti (2^]23|]. Extending the PIQMC method in this 
direction is left for future research. 



TABLE II. 


Path-integral Quantum Monte Carlo results for the , 


ground state energy of the hydrogen Hamiltonian, in units 


of q 2 /a (ao ■ 


= h 2 /Mq 2 ). The exact value is E = -0.5. 
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m 


Em 


20 


400 


-0.496 (± 0.004) 


20 


800 


-0.503 (± 0.005) 


40 


800 


-0.498 (± 0.006) 



TABLE III. Path-integral Quantum Monte Carlo results for the ground state energy of the helium Hamiltonian, in units of 
q 2 /ao- The experimental value is E = —2.904. 






m 




10 


400 


-2.84 (± 0.02) 


10 


800 


-2.88 (± 0.02) 


10 


1200 


-2.92 (± 0.03) 
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APPENDIX 



In PIQMC the simplest mehod for sampling paths is to change one degree of freedom at each Monte Carlo step. 
Usually this is rather inefficient and one adds Monte Carlo moves that make global changes of the path, e.g. moves that 
resembles the classical motion. In this appendix we present a more sophisticated scheme which we found performed 
very well at very low temperature. The basic idea is to change variables such that the kinetic energy term in the path 
integral becomes a diagonal quadratic form, i.e. 

m m 

J2(x k -x k+1 ) 2 = J2vl (28) 

fc=l k=2 

where x m+ i — x\. After some straightforward algebra one finds that the transformation from the {xj} to the {yj} is 
given by 

2 m-k+2 ( (to - k + l)x k -i + x m+1 \ 2 



m — k + 1 \ m— k + 2 

The expression for x k in terms of the {uj} reads 

Em — fc + l/m— j + l\ 1 ^ 2 , 
r—r[ -7—:) yj, 1 < k <m, 30 
. m-j + 1 \m-j + 2 J 

with x\ = y\. From (30) we conclude that the computational work for making a global change of the path (i.e. 
simultaneously changing all yj) is linear in m, hence optimal. It is also clear that the variable yi plays the role of the 
"classical" position. The variables j/2, • • ■ , Vm describe the quantum fluctuations. 



[1] K.E. Schmidt and D.M. Ceperley, in: The Monte Carlo Method in Condensed Matter Physics, ed. K. Binder, Springer, 
Berlin, 203 (1992). 

[2] H. De Raedt and W. von der Linden, in: The Monte Carlo Method in Condensed Matter Physics, ed. K. Binder, Springer, 

Berlin, 249 (1992). 
[3] H. De Raedt and M. Frick, Phys. Rep. 231, 107 (1993). 

[4] H. De Raedt, W. Fettes, and K. Michielsen, in: "Quantum Monte Carlo Methods in Physics and Chemistry", eds. M.P. 

Nightingale and C.J. Umrigar, NATO-ASI Series, 37 (Kluwer, The Netherlands 1999). 
[5] S. Lie and F. Engel, Theorie der Transformationgruppen, Teubner, Leipzig, 1888. 
[6] M. Suzuki, S. Miyashita, and A. Kuroda, Prog. Theor. Phys. 58, 1377 (1977). 
[7] H. De Raedt, and A. Lagendijk, Phys. Rep. 127, 233 7 (1985). 
[8] M. Suzuki, Phys. Lett. A201, 425 (1995). 

[9] R.P. Feynman and A.R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965. 
[10] H. Kleinert, Path integrals in Quantum Mechanics, Statistics and Polymer Physics, World Scientific, London, 1990. 
[11] I.H. Duru and H. Kleinert, Phys. Lett. B84, 30 (1979) 
[12] K. Symanzik, J. Math. Phys. 6, 1155 (1965). 
[13] H. De Raedt and B. De Raedt, Phys. Rev. A 28, 3575 (1983). 
[14] M. Suzuki, J. Math. Phys. 26, 601 (1985). 
[15] This is the case for all V(x) that are polynomial in x. 
[16] L.I. Schiff, Quantum Mechanics McGraw-Hill, New York, (1968). 
[17] G. Baym, Lectures on Quantum Mechanics, WA. Benjamin, Reading MA, (1969). 
[18] R.P. Feynman and H. Kleinert, Phys. Rev. A34, 5080 (1986). 
[19] H. Kleinert, Phys. Lett. B181, 324 (1986). 
[20] H. Kleinert, Phys. Lett. A118, 195 (1986). 
[21] W. Janke and B.K. Chang, Phys. Lett. B129, 140 (1988). 
[22] R. Giachetti and V. Tognetti, Phys. Rev. Lett. 55, 912 (1985). 
[23] R. Giachetti, V. Tognetti, and R. Vaia, Physica Scripta40, 451 (1989). 



7 



